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Abstract 

We study the critical behavior of frustrated systems by means of Pade-Borel 
resummed three-loop renormalization-group expansions and numerical Monte 
Carlo simulations. Amazingly, for six-component spins where the transition 
is second order, both approaches disagree. This unusual situation is analyzed 
both from the point of view of the convergence of the resummed series and 
from the possible relevance of non perturbative effects. 

P.A.C.S. numbers:05.50+q,05.70.Fh, 75.10.Hk,64.60.Cn, 75.10.-b 

Frustrated spin systems have been very much studied in their classical as well as quan- 
tum aspects. In particular, the critical behavior of 3D Stacked Triangular Antiferromagnets 
(STA) has deserved much attention |T|-|S| since, firstly it has many physical realizations in 
rare earth materials, secondly it is an archetype for frustrated systems, and thirdly it is 
directly related to the behavior of its two dimensional zero temperature, quantum counter- 
part. The frustration in this system comes from the fact that - for N > 1 component spin 
systems - the ground state is non collinear and shows the famous 120° structure. It is thus 
natural to believe that if the transition is second order, it belongs to a new universality class. 
Our present understanding of these systems comes as usual from theoretical renormalization 
group (RG) calculations, from Monte Carlo simulations and from experiments. The most 
impressive fact concerning these systems is that more than twenty years after the first works 
devoted to their study, there is still no agreement between these approaches. For instance, 
a calculation made in D = 4 - e predicts no stable fixed point for N in the interval: 
N c2 = 2.202 - 0.569e + 0.989e 2 < N < N c3 = 21.80 - 23.43e + 7.088e 2 and another made in 



D = 2 + e predicts a fixed point for any N larger than 2. Some experiments find a second 
order phase transition while others a weak first order. Moreover, the different approaches 
finding a second order transition do not find the same exponents, a fact that suggests that the 
theoretical or numerical approaches may miss a fundamental point (presence of topological 
defects, breakdown of perturbation theory, lattices of too small sizes. . . ). 

Our aim in this Letter is to shed light on this problem. We rely on the fact that the 
three-loop RG calculations made in D = 4 — e with e = 1 and directly in D = 3 find a critical 
value N C (D = 3) above which the transition is second order equal to 3.39 || and 3.91 [|J, 
respectively. Below it, it is supposed to be first order. We therefore expect a very weakly 
first order transition for N = 3 (and even perhaps for iV = 4) - a situation very difficult 
to test numerically. Therefore, instead of studying directly the physical, i.e. N = 3, spin 
system, we choose to study the following question: is there consensus between the results 
given by the usual RG approach based on the Landau- Wilson (0 4 - like) model and those 
obtained by Monte Carlo simulations for the values of D and N where a fixed point is found? 

Let us emphasize that the reliability of the RG approach for predicting the three- 
dimensional critical behaviors is not generic but has been demonstrated for most popular 
and simplest universality classes such as the usual 0(N) and cubic ones. The discrepancy 
between the perturbative results around d = 2 and d = 4 is in fact common to a wide 
class of systems among which the dipole locked phase of 3 He, electroweak phase transition, 
smectic liquid crystal, etc. Our study is therefore likely to be relevant to a much wider class 
of systems than the frustrated magnets. 

To tackle with our question, we study in d = 3 the largest possible N compatible with 
numerical possibilities where the usual recipes should work since in this case we are far 
above the line N C (D), the proximity of which could be the root of all the problems. Being in 
principle in the second order region, we expect to compute accurately the critical exponents 
both numerically and in a RG approach using the usual resummation techniques. The 
comparison between the results obtained by these two methods should be a test of the most 
powerful theoretical approaches in this non-ferromagnetic case. We also choose the value of 
TV such that the corresponding system does not show topological defects in order to eliminate 
a possible reason for the breakdown of perturbation theories. It turns out that iV = 6 is the 
ideal candidate. We present in this article numerical results for N = 6 as well as analytical 
ones for many N and compare them. 

Renormalization Group calculations. Let us first show our results obtained from the 
renormalization group analysis. The Landau- Wilson Hamiltonian relevant to the STA sys- 
tem is: 



H = - / d\r 



u j. ± * j. ± * , ™0 

2 



(1) 



The domain of parameters of interest is uq > and wq > 0. The calculations are based on 
the three-loop RG equations obtained earlier for the more complicated model having three 
independent quartic coupling constants @]. They are carried out directly in d — 3 and a 
Pade-Borel resummation of the relevant expansions is performed. Pade approximants [3/1] 
and [2/1] are employed for analytical continuation of the Borel transform series for the (3- 
functions and critical exponent 7, respectively. The Fisher exponent is evaluated by direct 
substitution of the fixed point coordinates into the corresponding expansion. 
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The fixed point of the RG flow diagram which controls the non-trivial (chiral) critical 
behavior is found. For N > 7 it turns out to be a stable node, for N = 5, 6, 7 this point 
manifests itself as a stable focus. The latter scenario looks quite new, i. e. is observed for the 
first time in STA systems, while the former one has been already discussed (see, e. g., Q). 
The estimates of critical exponents 7 and 77 for various iV are obtained from corresponding 
RG series. Making use of well-known scaling relations yields numerical values of the others. 
The results of our RG calculations are collected in Table 1 and presented, along with the 
other data, in Fig. 3. As is seen, critical exponents as functions of N demonstrate a cusp 
between N = 7 and N = 8 that reflects the abovementioned change of type of the fixed 
point governing the critical behavior. 

Monte Carlo results. We present now our MC simulations. We use six-component spins 
interacting via the Hamiltonian 

H = JijSi.Sj , (2) 

to) 

where the sum runs over all neighbors of the stacked triangular lattice (STA) and the 
interaction is chosen antiferromagnetic (J > 0). In the ground state the spins are planar with 
the three spins at the corners of each triangle forming a 120° structure. We use the standard 



Metropolis algorithm in combination with the over- relaxation algorithm |L0j . Between each 
Metropolis we use one over-relaxation step. This allows us to reduce the correlation time 
and obtain better statistics. For each size we use some hundred thousand steps to equilibrate 
our system and up to five millions steps to thermalize for the bigger sizes. We have repeated 
these simulations for different initial configurations (ordered or random) to make sure that 
our results do not depend on them. We use the histogram MC technique developed by 



Ferrenberg and Swendsen fLl| . From a simulation done at T , this technique allows us to 
obtain thermodynamic quantities at T close to T . We have studied our system in the finite 
size scaling (FSS) region |12|] and our simulations have been done at T s = 0.463. We consider 
L 2 * (2L/3) systems, where (L) 2 is the size of the planes, and 2L/3 is the number of planes 
in the z axis. First, to find the critical temperature T c , we use Binder's cumulant defined as 

U = 1- < M 4 > /(3 < M 2 > 2 ) (3) 

where the order parameter M is calculated in partitioning our lattice in three sublattices 
with only collinear spins and by summing each magnetization. We record the variation of U 
with T for various system sizes and then locate the intersection of these curves. We compare 
the values of U for two different lattice sizes L and L' = bL, making use of the condition 




U, 



bL 



u, 



1. (4) 



T=T n 



In Fig. [I], U is plotted as a function of the temperature for different sizes from L = 12 up to 
L = 36. Due to the presence of residual corrections to finite size scaling, one actually needs to 
extrapolate the results of this method for (In b)~ x — > 0. From these data, we extrapolate the 
value of T c (not shown) and obtain: T c = 0.4636(2) . We estimate the universal quantity U 
at T c (U*) as U* = 0.6545(15). With the value of T c we calculate the critical exponents using 
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log-log fit 0,0. We obtain from V x =< ME > / < M > - < E >, V 2 =< M 2 E > / < 
M 2 > — < E >, (Fig. g), estimates of 1/v, from the suceptibility \ = N < M 2 > /{k B T) 
(not shown) estimates of 7/z/, and from < M > (not shown) an estimate of (3 /v. Combining 
these results we obtain respectively, v = 0.700(11), 7/V = 1.975(20), and (3/u = 0.513(12). 
All our errors are calculated with the help of the Jackknife procedure |15[ and include the 
influence of the uncertainty in estimating T c . The results are summarized in Table [H] where 
the value of r\ has been calculated with the hyperscaling relation 7/V = 2 — rj. We note that, 
contrary to spins with two or three components, 77 is positive. This is due to the fact that 
for N = 6 the Renormalization Group flow is attracted by a true stable fixed point and not 
simply by a local minimum 

Discussion. The predictions of the renormalization-group ^-expansion technique for N = 
6-component spins listed in Table | do not agree with the Monte Carlo results given in Table 
|l[ We have plotted in figure |3] the results for v from the MC data, the RG ^-expansion 
(this work), the Local Potential Approximation method (LPA) [[]] and the 1/N expansion 
|T6| (first order). The six-loop RG results for the ferromagnetic case |T7| are also plotted for 



comparison. The interesting point is that these results are obtained by methods representing 
the state of the art in the field of critical phenomena. When applied to the systems belonging 
to the O(N) universality classes they indeed fit very well together. Let us emphasize that 
our numerical results are well converged and it seems unlikely that a non trivial behavior 
shows up at much larger system sizes THAT would resolve the discrepancy with the RG 
^-expansion results. 

In order to clear up what can be an origin of the marked discrepancy, we analyze the 
structure of the RG expansions employed. The main attention is paid to the vicinity of the 
chiral fixed point of the RG flow for N = 5, 6, 7 when this point is a focus. Contrary to 
the (unstable) fixed point governing the 0(./V)-symmetric critical behavior, the chiral point 
lies very close to the w axis being far from the u axis. In particular, for the case of interest 
N = 6 its coordinates are: u* = 0.0665, if* = 1.6025. In this region, the structure of the 
series of the /3-functions turns out to be unexpectedly irregular. As an example, we present 
here two "cuts" of the Borel-transformed expansion for /3 u (u, w) running through the chiral 
fixed point which clearly demonstrate such irregularity: 

P*{u, 1.6025) = -0.3607 + 0.7774m - 0.5004m 2 + 0.0339m 3 - 0.0055m 4 , (5) 
#f (0.0665, w) = 0.0643 - 0.0132m; - 0.1960m; 2 + 0.0346m; 3 + 0.0010m; 4 . (6) 

The coefficients in these formulas do not decrease monotonically with increasing their num- 
bers, and the expansion in powers of w is not alternating having coefficients with irregular 
signs. Therefore, the RG series for /3-functions would not demonstrate a good summability 
near the chiral fixed point and, subject to Pade-approximant-based analytical continuation 
and subsequent Borel integration, they are hardly believed to yield precise numerical results. 
Moreover, the Pade-Borel approximant for (3 U , taken at the chiral fixed point, as a function 
of the Borel variable t has a pole at t = 61.8 which is not dangerous practically but reflects 
the series poor summability. The difference between numerical results obtained within RG 
and MC approaches may be caused by an unfavorable structure of the RG expansions. On 
the other hand, for all N studied the chiral fixed point coordinate u* given by the resummed 
three-loop series remains positive preventing the RG expansions from losing Borel summa- 
bility in the domain of interest. Hence, fortunately, we do not face here this serious problem 
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that may spoil a perturbative RG analysis, as it occurs, say, when systems with quenched 



disorder are investigated fll8fl . This keeps calculations of the higher-order contributions to 
the RG functions of the model Eq.(JI|) meaningful and desirable. 

Can an account for four-loop or higher-order terms in the RG expansions significantly 
improve the situation? In principle, yes. The point is that the exact coordinates of the chiral 
fixed point may differ substantially from those given by the three-loop approximation and 
lie in the domain of the RG flow diagram where the perturbative expansions of /^-functions 
can be properly resummed. Higher-order terms added to the three-loop series may shift 
calculated fixed point coordinates toward their exact values thus making counterparts of the 
series (|5|-§|) better summable. To clear up whether such a situation really takes place, higher- 
order RG calculations have to be performed. Until this is done, an alternative scenario, i. 
e., the case when higher-order terms do not improve the summability, cannot be excluded. 

There is up to now only one other theoretical apparoach that allows quantitative calcu- 
lations in D = 3: the LPA method that consists in a truncation of the Wilson RG equations. 
Note that even if the LPA is missing the field renormalization and thus the anomalous dimen- 
sion rj, it is non perturbative since it is not based on a weak coupling expansion. However, 
although in our case the results obtained within this method are closer to the MC data than 
those obtained at the three-loop RG approximation, they show an unexpected dependence 
of v with N at small N. Moreover, used around d = 2, this approach contradicts the pertur- 
bative results obtained from the Non Linear Sigma model that are, in this dimension, well 
confirmed by simulations |20|| . They are anyway not enough accurate to draw a conclusion 
in d = 3. Since the LPA is known to be the first order of a systematic derivative expansion 
of the effective action, it is desirable that the next order be computed. 

Let us now remark that even if the three-dimensional physics was well reproduced by 
our calculations, it would remain that a coherent picture of the critical thermodynamics of 
the frustrated systems would require to understand the discrepancy between the NLcr model 
approach and the Landau- Wilson one. A striking difference between both approaches is that 
near two dimensions the low temperature expansion of the NLcr model predicts that a new 
"current" term of the form (c/>*Vc/>) 2 is relevant 0. This term appears to be fundamental 
since for N = 3 it allows to find a fixed point with an 0(4) symmetry. Being highly non 
renormalizable near four dimensions, it is irrelevant and forgotten. There is thus another 
scenario than the numerical unreliability of the three-loop RG approximation, namely, that 
the Landau- Wilson Hamiltonian Eq.fll]) itself is incomplete in three dimensions (remember 
that the presence of topological defects cannot be invoked here since there are no such 
defects for N = 6). As it was suggested for the Abelian Higgs transition, this could be 
interpreted as the necessity to have recourse to the NLcr model description and to abandon 
that of the Landau- Wilson model. Note, however, that it is very doubtful that the analysis 
made around two dimensions can be extended straightforwardly for any N up to d = 3 
since i) for iV = 3-component spins the 0(4) fixed point found in D = 2 + e has been 
shown to disappear in a non trivial dimension strictly smaller than three in a closely related 
model - the principal chiral model |HJ - and since ii) an 0(4) behavior has neither been 
seen experimentally nor numerically for N = 3 and d = 3. Thus, the perturbative analysis 
of the NLcr model fails also above two dimensions. However, it remains that a coherent 
picture of the behavior of frustrated systems for all N and d should include the results of 
the NLcr model and therefore explain why and when the current term starts to be relevant 
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as a function of N and d. If this happens to be around d = 3 for N ~ 0(1), it could perturb 
the RG (/-expansion results presented here and explain why otherwise powerful methods do 
not work properly in our case. In any case, we believe that our results for N = 6 constitute a 
clear challenge for the theoretical approaches which is perhaps not out of reach from higher 
order calculations and/or improvement of the LPA method. 
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TABLE I. Critical exponents calculated by RG 
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TABLE II. Critical exponents obtained by MC. Calculated with a = 2 — dv. Calculated with 
r] = 2 — j/u. 
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FIGURES 
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FIG. 1. Binder's parameter U as function of the temperature for different sizes L (in the left 
part of the figure from down to up L=12, 15, 18, 21, 24, 27, 30, 36). The arrows show the estimated 
critical temperature T c and the temperature of our simulations T s . 




FIG. 2. Values of V\ and V2 as function of L in a ln-ln scale at T c . The value of the slopes 
gives \jv and we obtain v= 0.698(12) for V x and 0.702(13) for V 2 . The smallest size (L = 12) is 
not included in our fits. 
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FIG. 3. v for different methods (see text). 
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